home *** CD-ROM | disk | FTP | other *** search
/ MacWorld 1997 September / Macworld (1997-09).dmg / Serious Software / Cherwell Scientific Demos / pro Fit / pro Fit 5.0 demo (fpu).sea / pro Fit 5.0 demo (fpu) / Functions & Programs / Plot Vector Field < prev    next >
Text File  |  1996-04-15  |  3KB  |  96 lines
  1. {
  2. This program plots a vector field. The input data must be arranged in
  3. four columns:
  4.  the x- and y-coordinates of the vectors (two columns)
  5.  the length of each vector
  6.  the angle of each vector (in radian)
  7. }
  8.  
  9. program VectorField;
  10.  
  11. const 
  12.       arrowAngle = 30;{ arrow opening angle }
  13.  
  14. var 
  15.      sinA, cosA: extended; { sine and cos of arrow angle times arrow length }
  16.            xmin,xmax,ymin,ymax:extended;    { ranges for the new graph }
  17.            arrowLength:extended; 
  18.            xCol,yCol,magnitudeCol,AngleCol:integer; { data columns }
  19.      i:integer;
  20.  
  21. procedure initialize;
  22. begin
  23.     xCol:=1;yCol:=2;magnitudeCol:=3;angleCol:=4;
  24.  arrowLength:=0.5;
  25. end;
  26.  
  27. procedure SetRanges;
  28. var newxMin,newyMin,xVal,yVal:extended;
  29. begin
  30.  xMin:=inf;xMax:=-inf;yMin:=inf;yMax:=-inf;
  31.  for i:=1 to NrRows do
  32.   if DataOK(i,xCol) and DataOK(i,yCol)and DataOK(i,magnitudeCol) and DataOK(i,AngleCol) then
  33.   begin
  34.    xVal:=data[i,xCol];yVal:=data[i,yCol];
  35.    if xVal<xMin then xMin:=xVal;
  36.    if yVal<yMin then yMin:=yVal;
  37.    if xVal>xMax then xMax:=xVal;
  38.    if yVal>yMax then yMax:=yVal;
  39.   end;
  40.   newxMin:=xMin-(xMax-xMin)/10;
  41.   xMax:=xMax+(xMax-xMin)/10;
  42.   xMin:=newxMin;
  43.   newyMin:=yMin-(yMax-yMin)/10;
  44.   yMax:=yMax+(yMax-yMin)/10;
  45.   yMin:=newyMin;
  46. end;
  47.  
  48.  
  49. procedure DrawArrow(x,y, vx,vy);
  50.   { draws an arrow at position x,y with vector vx, vy }
  51.   var tipX, tipY: real;
  52.       ax, ay, r: real;
  53. begin
  54.   r := sqrt(sqr(vx) + sqr(vy));
  55.   if r > 1e-30 then
  56.   begin
  57.    tipX := x + vx/2;
  58.    tipY := y + vy/2;
  59.    moveto(x-vx/2, y-vy/2);
  60.    lineto(tipX, tipY);
  61.    if r >= arrowLength then   { show arrow only if the shaft is long enough }
  62.    begin
  63.     vx := vx/r; vy := vy/r;
  64.     ax := cosA*vx - sinA*vy;
  65.     ay := sinA*vx + cosA*vy;
  66.     moveto(tipX-ax, tipY-ay);
  67.     lineto(tipX, tipY);
  68.     ax := cosA*vx + sinA*vy;
  69.     ay := -sinA*vx + cosA*vy;
  70.     lineto(tipX-ax, tipY-ay);
  71.    end;
  72.   end;
  73. end; { DrawArrow }
  74.  
  75. begin
  76.         Input('$Cx-coordinates:', xCol,
  77.               '$Cy-coordinates:', yCol,
  78.               '$Clengths:', magnitudeCol,
  79.               '$Cangles [rad]:', angleCol,
  80.               'length of arrows', arrowLength);
  81.   sinA := sin(arrowAngle*π/180) * arrowLength;
  82.      cosA := cos(arrowAngle*π/180) * arrowLength;
  83.         SetRanges;
  84.   CreateNewGraf(xMin,xMax,yMin,yMax,0,0); 
  85.  
  86.   OpenCurve('vector field');
  87.      for i:=1 to NrRows do
  88.       if DataOK(i,xCol) and DataOK(i,yCol)and DataOK(i,magnitudeCol) and DataOK(i,AngleCol) then
  89.       begin
  90.                 DrawArrow(data[i,xCol],data[i,yCol],
  91.                 data[i,magnitudeCol]*cos(data[i,AngleCol]),data[i,magnitudeCol]*sin(data[i,AngleCol]));
  92.             end;
  93. end;
  94.  
  95.   
  96.